Map with base graphics

Let’s create a very basic “map”, showing the location of acoustic samples along the West Coast. The first uses base graphics (yuck!) and a cartesian coordinate system.

# Load some ship track data
cps <- read_csv(here("Data/cps_nav.csv"))

# A basic 
plot(cps$long, cps$lat)

# plot(cps$long, cps$lat, cex = cps$cps.nasc)
# plot(cps$long, cps$lat, cex = sqrt(cps$cps.nasc/1000))

Map using ggplot2

Here are the same data plotted using ggplot2.

# A very basic ggplot2 scatter plot of the same data
ggplot(cps, aes(long, lat)) + geom_point()

ggmap

ggmap is a mashup of ggplot and online mapping services. There are numerous map sources available to ggmap (e.g., Google, OpenStreetMaps, Stamen). We’ll look quickly at Stamen and Google maps. Withing each of those, there are multiple options (e.g., Google has satellite, terrain, etc.). Each has different arguments to specify the map boundaries.

# Set west coast boundaries for stamen maps
wc.bounds.stamen <- c(left = min(wc.long), bottom = min(wc.lat),
                      right = max(wc.long), top = max(wc.lat))
# Download stamen map of west coast; zoom = 6 seems good
wc.map.stamen.toner <- get_stamenmap(wc.bounds.stamen, zoom = 6, maptype = "toner-lite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw()
# Display the stamen map
wc.map.stamen.toner

# Add layers to map
wc.stamen1 <- wc.map.stamen.toner + 
  geom_point(data = cps, aes(long, lat, group = transect, colour = factor(transect), size = cps.nasc),
             show.legend = F) +
  geom_text(data = locations, aes(long, lat, label = name)) +
  ggtitle("Basic options")

wc.stamen2 <- wc.map.stamen.toner + 
  geom_point(data = cps, aes(long, lat, group = transect, colour = factor(transect), size = cps.nasc),
             show.legend = F) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text(data = locations, aes(long, lat, label = name), 
            colour = "red", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) + 
  ggtitle("Formatted labels")

wc.stamen3 <- wc.map.stamen.toner + 
  geom_point(data = cps, aes(long, lat, group = transect, colour = factor(transect), size = cps.nasc),
             show.legend = F) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text_repel(data = locations, aes(long, lat, label = name), 
                  size = 2, segment.colour = "black", segment.alpha = 0.5) +
  ggtitle("Formatted and repelled lables")

# Combine maps
wc.grid.stamen <- plot_grid(wc.stamen1, wc.stamen2, wc.stamen3, nrow = 1)

# Save map
ggsave(wc.grid.stamen, filename = here("Figs/wc_map_stamen_toner.png"), width = 10, height = 7)
# Print map
include_graphics(here("Figs/wc_map_stamen_toner.png"))

# Reduce label list
label.list <- c("Monterey Bay","San Francisco","Cape Flattery","Crescent City",
                "Newport","Point Conception","Cape Mendocino","Columbia River",
                "Cape Blanco","Bodega Bay","Westport","Fort Bragg",
                "Morro Bay","Long Beach","Cape Scott","San Diego")

locations <- filter(locations, name %in% label.list)

# Set west coast boundaries for stamen maps
wc.bounds.google <- c(mean(wc.long),mean(wc.lat))
# Download stamen map of west coast
wc.map.google <- get_map(location = wc.bounds.google, zoom = 4, source = "google", maptype = "satellite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
  xlim(wc.long) + ylim(wc.lat)

wc.google1 <- wc.map.google +
  geom_point(data = cps, aes(long, lat, group = transect, colour = cps.nasc, size = cps.nasc)) +
  scale_colour_gradientn(colors = rev(rainbow(10))) +
  geom_point(data = locations, aes(long, lat), colour = "white") +
  geom_shadowtext(data = locations, aes(long, lat, label = name), 
                  colour = "white", bg.color = "black", size = 3, 
                  nudge_x = 0.25, hjust = 0, angle = 45, fontface = "bold.italic") +
  ggtitle("Shadow text, scale color & size")

wc.google2 <- wc.map.google +
  geom_point(data = cps, aes(long, lat), colour = "gray50", size = 0.25, alpha = 0.75) +
  geom_point(data = filter(cps, cps.nasc > 0), 
             aes(long, lat, group = transect, colour = cps.nasc, size = cps.nasc), alpha = 0.75) +
  scale_colour_viridis_c(option = "magma") +
  geom_point(data = locations, aes(long, lat), colour = "white") +
  # Configure legend guides
  guides(colour = guide_legend(), size = guide_legend()) +
  geom_shadowtext(data = locations, aes(long, lat, label = name), 
                  colour = "white", bg.color = "black", size = 2, 
                  nudge_x = 0.25, hjust = 0, angle = 45, fontface = "bold.italic") +
  ggtitle("Viridis colours, grouped legend")

# Arrange plots
wc.grid.google <- plot_grid(wc.google1, wc.google2, nrow = 1)

# Save map
ggsave(wc.grid.google, filename = here("Figs/wc_map_google.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_google.png"))

# Set west coast boundaries for stamen maps
mb.bounds.google <- c(mean(mb.long),mean(mb.lat))
# Download stamen map of west coast
mb.map.google <- get_map(location = mb.bounds.google, zoom = 10, 
                         source = "google", maptype = "satellite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
  xlim(mb.long) + ylim(mb.lat)

mb.google1 <- mb.map.google +
  geom_point(data = filter(locations, name != "Monterey Bay"), aes(long, lat), colour = "white") +
  geom_text(data = filter(locations, name != "Monterey Bay"), aes(long, lat, label = name), 
            colour = "white", size = 4, vjust = 0, nudge_y = 0.01)
# Save map
ggsave(filename = here("Figs/mb_map_google.png"), width = 7, height = 7)
# Print map
include_graphics(here("Figs/mb_map_google.png"))

rnaturalearth

rnaturalearth for shoreline and country files. Below is a comparison of the medium- and large-scale country datasets.

# Import NaturalEarth coastlines and countries
coast_sf        <- ne_coastline(scale = "medium", returnclass = "sf")
coast_sf_lg     <- ne_coastline(scale = "large", returnclass = "sf")
countries_sf    <- ne_countries(scale = "medium", returnclass = "sf")
countries_sf_lg <- ne_countries(scale = "large", returnclass = "sf")

# Display countries
plot(countries_sf)

# View unique regions
sort(unique(countries_sf$region_wb))

## [1] "Antarctica"                 "East Asia & Pacific"       
## [3] "Europe & Central Asia"      "Latin America & Caribbean" 
## [5] "Middle East & North Africa" "North America"             
## [7] "South Asia"                 "Sub-Saharan Africa"

sort(unique(countries_sf$subregion))

##  [1] "Antarctica"                "Australia and New Zealand"
##  [3] "Caribbean"                 "Central America"          
##  [5] "Central Asia"              "Eastern Africa"           
##  [7] "Eastern Asia"              "Eastern Europe"           
##  [9] "Melanesia"                 "Micronesia"               
## [11] "Middle Africa"             "Northern Africa"          
## [13] "Northern America"          "Northern Europe"          
## [15] "Polynesia"                 "Seven seas (open ocean)"  
## [17] "South-Eastern Asia"        "South America"            
## [19] "Southern Africa"           "Southern Asia"            
## [21] "Southern Europe"           "Western Africa"           
## [23] "Western Asia"              "Western Europe"

# Select only certain regions to speed plotting
na_sf    <- filter(countries_sf, subregion %in% c("Northern America","Central America")) 
na_sf_lg <- filter(countries_sf_lg, subregion %in% c("Northern America","Central America")) 

# Create West Coast map (medium scale)
wc.ne <- ggplot() + 
  geom_sf(data = na_sf) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text(data = locations, aes(long, lat, label = name), size = 4) +
  scale_x_continuous(name = "Longitude", limits = wc.long) + 
  scale_y_continuous(name = "Latitude", limits = wc.lat) +
  theme_bw() + 
  coord_sf(xlim = c(-123.5,-121.5), ylim = c(36.5,38.5)) +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5),
          legend.position =  c(0,0),
          legend.justification = c(0,0),
          legend.background = element_blank(),
          legend.key = element_blank(),
          panel.background = element_rect(fill = alpha("lightblue", 0.5)),
          plot.title = element_text(hjust = 0.5),
          panel.grid.major = element_line(color = "white")) +
  ggtitle("Medium scale") 

# Create West Coast map (large scale); and with shadowtext labels
wc.ne.lg <- ggplot() + 
  geom_sf(data = na_sf_lg) +
  geom_point(data = locations, aes(long, lat)) +
  geom_shadowtext(data = locations, aes(long, lat, label = name), colour = "white",
                  size = 4, bg.colour = "black") +
  scale_x_continuous(name = "Longitude", limits = wc.long) + 
  scale_y_continuous(name = "Latitude", limits = wc.lat) +
  theme_bw() + 
  coord_sf(xlim = c(-123.5,-121.5), ylim = c(36.5,38.5)) +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5),
          legend.position =  c(0,0),
          legend.justification = c(0,0),
          legend.background = element_blank(),
          legend.key = element_blank(),
          panel.background = element_rect(fill = alpha("lightblue", 0.5)),
          plot.title = element_text(hjust = 0.5),
          panel.grid.major = element_line(color = "white")) +
  ggtitle("Large scale")

# Create map grid
wc.grid.ne <- plot_grid(wc.ne, wc.ne.lg, nrow = 1)

# Save map
ggsave(wc.grid.ne, filename = here("Figs/wc_map_natEarth.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_natEarth.png"))

Merge

Buffer

Compute area

etc.

elevatr

The elevatr package provides easy access to elevation data from a variety of sources. Not fully functional at the moment…

FedData

FedData Github repo and tutorial with examples of various data sources (e.g., .

FedData Tutorial

Presently not executed; still under development

# Get a map of Moss Landing
ML.bbox <- get_map("Moss Landing, CA", zoom = 14) 
# Plot the map
ggmap(ML.bbox) 
# Get the bounding box from its attributes
bb <- attr(ML.bbox, "bb")
# Create an extent polygon from the Monterey Bay bounding box
ML.extent <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon,
                                            bb$ll.lat, bb$ur.lat), 
                             proj4string="+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
# Map the bounding box to check the results
ggmap(ML.bbox) + geom_polygon(aes(x=long, y=lat), size=3, color='purple',
                             data=ML.extent, fill=NA)
# Download National Elevation Data
ned_ML <- get_ned(template = ML.extent, label="ned_ML", 
                    res="1", force.redo = F, raw.dir = here("Data"))
# Convert raster to data frame for plotting in ggplot
ned_df <- rasterToPoints(ned_ML) %>% 
  data.frame() %>% 
  rename(lat = y,
         long = x,
         elev = grdn37w122_1)

# Map elevation
ggplot(ned_df, aes(long,lat,fill = elev)) + geom_raster() +
  geom_contour(aes(z = elev), alpha = 0.5, binwidth = 1) + 
  scale_fill_viridis_c(name = "Elevation (m)") +
  theme_bw()

Movement analysis

A collection of movement analysis packages